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Abstract 

Vocal fold (VF) motion is fundamental to voice production and di- 
agnosis in speech and health sciences. The motion is a consequence of 
air flow interacting with elastic vocal fold structures. Motivated by ex- 
isting lumped mass models and known flow properties, we propose to 
model the continuous shape of vocal fold in motion by the two dimen- 
sional compressible Navier-Stokes equations coupled with an elastic 
damped driven wave equation on the fold cover. In this paper, instead 
of pursuing a direct two dimensional numerical simulation, we derive 
reduced quasi-one-dimensional model equations by averaging two di- 
mensional solutions along the flow cross sections. We then analyze the 
oscillation modes of the linearized system about a flat fold, and found 
that the fold motion goes through a Hopf bifurcation into temporal os- 
cillation if the flow energy is sufficient to overcome the damping in the 
fold consistent with the early models. We also analyze the further re- 
duced system under the quasi-steady approximation and compare the 
resulting vocal fold equation in the small vibration regime with that 
of the Titze model. Our model shares several qualitative features with 
the Titze model yet differs in the specific form of energy input from 
the air flow to the fold. Numerical issues and results of the quasi-one- 
dimensional model system will be presented in part II (view resulting 
web VF animation at http: / /www.ma.utexas.edu / users /jxin| ) . 
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1 Introduction 



Vocal folds consist of two lips of opposing ligaments and muscle at the top 
of the trachea and join to the lower vocal tract. When air is expelled at suf- 
ficient velocity through this orifice (the glottis) from the lung after a critical 
lung pressure is achieved, the folds vibrate and act as oscillating valves in- 
terrupting the airflow into a series of pulses. These pulses of air flow serve as 
the excitation source for the vocal tract in all voiced sounds. For background 
description of vocal folds and speech production, see [f| , || , p3[ among oth- 
ers. Mathematical modeling of vocal fold motion can help us understand 
the voiced sound generation and make synthetic machine voice more natu- 
ral. Such an understanding also helps the medical diagnosis and correction 



of voice disorders ||, |T5| . 



Vocal fold modeling relies on our knowledge of aerodynamics and biome- 
chanics, as the problem is basically airflow through a flexible channel with 
time-varying cross section. One of the best known models of vocal folds is 
the two mass model by Ishizaka and Flanagan ||. The vocal fold is modeled 
by upper and lower masses (rrii, i = 1,2) connected by an elastic spring and 
attached to the wall by upper and lower springs with damping. The two 
mass motion is given by the ODE system: 

miX^u + + hxi + k c (-l) l+1 (x 1 - x 2 ) = F h i = 1,2, (1.1) 

where r^s are the damping coefficients; k^s, k c , are the elastic restoring con- 
stants (stiffness coefficients) for the upper, lower and side springs respectively; 
Xi s are the displacements of two masses from their prephonatory equilibrium 
positions; Fj's are driving forces acting on the two masses from glottal air 
pressure. To close the equations and complete the model, simplified as- 
sumptions are made to allow a calculation of Fi from the fluid flow so that 
Fi = Fi(xi,x 2 , di, k, P s ), where P s is the lung pressure, di the thickness of rrii 
along the flow direction, U the length of m; transverse to the flow direction. 
The flow is assumed to be one dimensional, inviscid, and quasisteady so that 
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Bernoulli's law can be applied. Moreover, F{ and elastic forces kiXi need to 
be modified empirically when the vocal folds tend to a closure to facilitate a 
pressure buildup to reopen the folds. In spite of these oversimplifications, the 
two mass model is able to capture many aspects of oscillation features. Small 
amplitude oscillation appears as a Hopf bifurcation from a steady state as 
P s passes a threshold value ||. More recently, Lucero [12| studied large am- 



plitude oscillations and found coexistence of multiple equilibria and studied 
their stability and bifurcations. 

The two mass model has been widely used as a glottal source for speech 
synthesis, and many improvements and extensions have been made to incor- 
porate additional physics. For example, see |§ for including effects of flow 
viscosity and flow separation; see |2(|, [26], for works on four and multiple 



mass models and applications. In M , a continuum elastic system under nat- 
ural boundary condition is used to find intrinsic vibrational eigenmodes of 
vocal fold tissues. 

Though two mass model is already simple looking, the theshold oscilla- 
tion condition on minimum lung pressure is still implicit analytically. Titze 
proposed his celebrated body-cover model in 1988 [^2| . The model is based 
on the hypothesis that during fold oscillation, the fold cover (epithelium and 
superficial layers of vocal ligament) propagates a surface mucosal wave in 
the direction of the airflow, and the body (deep layer of the vocal ligament 
and muscle) is stationary. The fold shape is approximated as a straight line 
connecting the fold entry of height hi and the fold exit of height h.2- Taylor 
expanding a mucosal wave with constant velocity c, Titze |22| approximated: 



h = 2(h + h + rht), h 2 = 2(h + h- rh t ), (1.2) 

where r = 2L/c the travel time for the mucosal wave to reach exit from en- 
trance, L half length of glottis along the flow direction, ho half glottal width, 
h fold displacement. The fold motion is lumped onto the fold midpoint, and 
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postulated as: 

1 r L 



1 

mh u + rh t + kh = — / P(x)dx = P g , (1.3) 

ALj J —L 

where m, r, and k are the mass, damping, and stiffness of the oscillating 
portion of the vocal fold about the midpoint, P(x) is the pressure distribution 
along the flow direction x. Using Bernoulli's law and linear fold shape, pZ 
showed in particular that for small fold vibration about rectangular fold 
equilibrium position h : 

P ' = w-$ = T7ITT- (L4) 

«i h + h + Tht 

where P s is the subglottal pressure, and exit pressure is zero. Linearizing 
( |1.4|) about h , one has: 

mh tt + (r - 2rP s h 1 )h t + kh = 0, (1.5) 

from which follows the threshold pressure condition setting the effective 
damping coefficient to zero: 

P s> * = h r/2r. (1.6) 

So the oscillation threshold pressure is proportional to damping constant r, 
and glottal half width /i - 

Titze model explains in a transparent fashion the formation of oscillation, 
and provides a very simple onset condition (|1.6|) . It shows how the energy 
input from the airflow balances the intrinsic fold damping due to the h t 
dependence of the fold driving force P g . The agreement of Titze model and 
the two mass model is discussed in |p2| , f23f . The onset condition (|L 



especially P s ^ being linear in r, is also supported by experiments P3| |. For 
results of Titze model on converging and diverging prephonatory fold shapes, 
large amplitude oscillations, and hysteresis phenomenon at onset-offset, see 

a, 0, h. 
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In either the two mass model or Titze model, the vocal fold shape is not 
described as a continuous curve as we see in physical reality, and that the 
treatment of air flow is rather crude. It is our goal to seek a more accurate 
and more first principle based model so that these two aspects are much 
improved. The fluid characteristics of vocal flow during phonation have been 
noted in ||, [^2j, [p3[| , [ |i~2|l , [|i~4}1 , among others, and also recently studied 



by Pelorson et al [T7[] using unsteady-flow measurements and visualization 
techniques. The flow is essentially two dimensional with Mach number on 
the order of 10 _1 , and Reynolds number on the order of 10 3 . The flow in the 
bulk is nearly inviscid except in the viscous boundary layers. 

Our starting point is the two dimensional compressible isentropic Navier- 



Stokes equations for the air flow. Motivated by |8j and ||22|| , we model the 
vocal fold as a finite mass elastic tube of cross sectional area A(x, t) with 
elastic attachments (muscles) onto nearby walls (bones). We shall either 
take rectangular cross section (length 2L in x, width 2w in z, and variable 
height 2h in y, then A = 4wh) or an elliptical shaped cross section (with 
principal axes of lengths 2h in y, and 2w in z). The air flows from x = — L 
to x = L, symmetric across y = 0, and is independent of z. The purpose 
of keeping the z dimension is for later reconstruction of a three dimensional 
picture of fold motion in part II @. Assuming that the flow is essentially 
along x direction, and its variation transverse to the tube is small except in 
the boundary layers near the folds, we average the flow on each cross section 
and derive a reduced quasi-one-dimensional system. The motion of fold cover 
is modeled by a damped driven elastic wave equation, considering the elastic 
tension on the fold cover and elastic forces in the attached muscles. 

The reduced quasi-one- dimensional aerodynamic equations are: 
• conservation of mass: 

{Ap) t + {puA) x = 0, (1.7) 

p air density, u air velocity; 
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• conservation of momentum: 

u t + uu x = — p x + + ^A-\Au x ) x - -fA- 1 A tx , 1.8 
p A 3p 3p 

p air pressure, p air viscosity. 

The force balance on the fold cover gives: 

• the dynamic boundary motion: 

m{A - A eq ) tt = <x(A - A eq ) xx - a{A - A eq ) t - (3{A - A eq ) + Sp + f m . (1.9) 

Here: m is the fold mass density; a is the longitudinal elastic tension of the 
fold ||, a is the muscle damping coefficient; (3 is an elastic modulus 
modeling the vibration property of the fold in the transverse; S is a cross 
section shape factor, S = Aw for rectangular cross sections, and S = nw for 
elliptical cross sections; f m is a prescibed function to model muscle tone so 
that a particular fold shape A eq (converging or diverging or flat) serves as 
an equilibrium state. In general, a, f3, a are functions of x, to model the 
varying stiffness of fold, since the fold cover is stiffer towards x = ±L. 
We shall write (11.91) into: 



mA a = oA xx - ocA t - (3 A + Sp + f m , (1.10) 

where f m = f m (x) is prescribed forcing. 
• The equation of state: 

p = Kp^, 7 > 1, K > 0. (1-11) 

The system (|1.7| )- (|1.11|) is closed, and we solve an initial boundary value 
problem on x G [— L, L] with proper in- flow/out-flow boundary conditions 
and initial fold shape. We shall consider the inviscid limit of our model by 
setting p = and normalize m — 1 while analyzing the oscillation onset in 
this paper. We shall show via a stability analysis that the linearized system 
near a flat fold admits (fold cover) waves in A of the form sin(kx — cut), for k 
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and uj real and nonzero if the glottal airflow has enough energy to overcome 
fold damping. This is the kind of fold surface mucosal wave hypothesized 
in Titze model [ 2^] . We also obtain similar onset conditions for small am- 
plitude oscillations due to the occurence of a pair of imaginary eigenvalues, 
consistent with early models. Under the quasi-steady approximation as in 



Titze model, we obtain a single fold equation resembling (1.5), the correction 
term is nonlocal but does depend on ht, and it plays the role of transferring 
aerodynamic energy onto the fold. 

We remark that the viscous effect is important when the fold is near 
closure, or when the fold is diverging enough so that flow separation occurs 



inside the glottis f22|, [[17]], ||. Our model above has not yet taken this 
aspect into full consideration. However, since flow separation decreases the 
pressure (pressure is zero from the point where flow detaches from the fold ) 
towards the exit x = L, and the fold is relatively stiff there, we can minimize 
the effect on fold motion by modeling the increasing stiffness of fold near 
the exit. Our model captures the compressibility of air flow which is critical 
during the opening of the folds jl6| . The viscous effect during the opening 
or near closure of vocal fold is delicate and requires further modeling, yet we 
have not found the lack of it to be a major problem in our simulation of the 
opening and closing of fold motion in part II [[5J. 

In studying collapsible tubes (|y|, [ll|]) and air flow through duct of 
spatially varying cross section fl25 |, it is common to use ( |1.7| ) with the one 
dimensional unsteady Euler equation which is ( |l.<j| ) with the last three terms 
omitted. The major differences in the two modeling problems are: (1) a 
vocal fold is fast oscillatory in time (e.g. 100 - 200 Hz), (2) the vocal fold 
carries mass, and a dynamic (damped driven wave) equation is necessary to 
describe the fold motion; moreover, the vocal fold has mechanical damping. 
In contrast, collapsible tubes are massless and damping free JT^j. The con- 
servation of mass and momentum in the collapsible tube system is recovered 
when A t <C A and \i — > 0. Though in both collapsible tube and vocal flow 
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problems, the flow is close to being incompressible, we opt not to make such 



an approximation as in fll9| . This is because an equation of state has to be 



determined in lieu of the natural equation of state ( |1.11|) when we consider 
quasi-one-dimensional reduction. For collapsible tubes, see |19[] and |TT] , the 
tube cross section is related to the pressure p by a tube law: p = A ni — A n2 
with p = constant. The tube law needs experimental measurements and de- 
pends also on materials involved. Due to lack of our knowledge of whether 
such a law exists for vocal fold, we choose to use the natural equation of state 
for air ( |1 . 1 1|) and work with compressible flow equations. It is interesting to 
investigate how good it is to use the incompressible Navier-Stokes equations 
for the two dimensional flow as a viable alternative. To rephrase our formu- 
lation, the p and p are related by the equation of state, the gamma gas law 
fll.ll ); then the cross section A is related to p dynamically. 

The rest of the paper is organized as follows. In section 2, we derive 
the fluid part of the quasi-one- dimensional system (L7|)- ([L.11|) from the two 
dimensional isentropic compressible Navier-Stokes equations by averaging 
solutions over the flow cross section. In section 3, we perform a linear stability 
analysis near a flat (rectangular) fold and deduce onset condition for the 
minimum flow energy to ensure the occurence of a pair of pure imaginary 
eigenmodes. In section 4, we make the quasi-steady approximation within 
our model, derive the single equation for the fold cover, obtain an onset 



condition for small amplitude oscillation, and compare with findings in |22 
Numerical issues, and numerical results will be presented in part II M. 



2 Derivation of Reduced Flow Equations 

We derive the fluid part of the model system assuming that the fold varies 
in space and time as A = A(x,t). Consider a two dimensional slightly vis- 
cous subsonic air flow in a channel with spatially temporally varying cross 
section in two space dimensions, Q = fi (t) = {(x,y) : x G [—L,L],y G 
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[—A(x, t)/2, A(x, t)/2]}, where A(x, t) denotes the channel width with a slight 
abuse of notation, or cross sectional area since the third dimension is uni- 
form. The two dimensional Navier-Stokes equations in differential form are 
(0, page 147): 

• conservation of mass: 

p t + V -{pu) = 0; (2.1) 

• conservation of momentum: 

(pu)t = -V- (p(u® u)) + div(a); (2.2) 
where a is the stress tensor, a = (a^) = —pSij + d^, and: 

dij = 2/i(ey - -^P<%), <;., :y (''/..,-, • Hj.,A- {xi,x 2 ) = {x,y); 

H is the fluid viscosity; Q(t) is any volume element of the form: 

n(t) = {( Xl y):xe [a, b] C [-L, L],y G [-A(x, t)/2, A(x, t)/2}.}. (2.3) 

The equation of state is Ql.llj) . 

The boundary conditions on (p, u) are: 

(1) on the upper and lower boundaries y = ±A(x,t)/2, p y = 0, and u = 
(0, ±A t /2), the velocity no slip boundary condition; 

(2) at the inlet, u(—L,y,t) = ui, a prescribed inlet velocity, p(—L,y,t) = pi, 
a prescribed inlet density (deduced from input pressure); 

(3) at the outlet, u x (L,y,t) = 0. 

We are only concerned with flows that are symmetric in the vertical. 
For positive but small viscosity, the flows are laminar in the interior of Qq 
and form viscous boundary layers near the upper and lower edges. The 
vertically averaged flow quantities are expected to be much less influenced 
by the boundary layer behavior as long as A(x, t) is much larger than 0(p 1 ^ 2 ). 
We also ignore effects of possible flow seperation inside Qq when it becomes 
divergent with large enough opening. 
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Let us assume that the flow variables obey: 



\ u i,y\ ^ \ u i,x\, \ u 2, y \ "C |ui,a;|) away from boundaries of fl , 
\u y \ \u x \, near the boundaries of Qq, 

\Py \ ^ I Pa; I) throughout fi . (2-4) 

These are consistent with physical observations in the viscous boundary lay- 
ers (Hj]], page 302), namely, there are large vertical velocity gradients, yet 
small pressure or density gradients in the boundary layers. The boundary 
layers are of width 0(p 1 ^ 2 ). Denote by p, Hi, the vertical averages of p and 
Ui. Note that the exterior normal n = (—A x /2, 1)/(1 + A 2 X /A) 1/2 if y = A/2, 
ft = (~A x /2, -1)/(1 + Al/Af' 2 if y = -A/2. 
Let a = x, b = x + 5x, 5x <C 1, t slightly larger than to- We have: 

4- / p dV = ^[ P J(t)dV=( Pt J(t)dV+ [ pJ t dV, (2.5) 
at Jn(t) at yn(to) y«(to) ^o(to) 



where J(t) is the Jacobian of volume change from a reference time to to t 

A(t ) 



Since f2(t) is now a thin slice, J(t) = 4rpr for small <5x, and J t = A t (t)/A(to). 



The second integral in (|2.5|) is: 



/ P J t dV = p^-A(t )6x = pA t (t)6x. (2.6) 



The first integral is simplified using (p.l[) as: 



Pt J(t)dV = Pt dV=- pu-ndS. (2.7) 

a(t ) jn(t) ^9n(t) 

We calculate the last integral of ( |2.7| ) further as follows: 

rA/2 [A/2 

pu-nds = / (—pui)(x,y,t)dy+ / (pui)(x + 5x,y,t) dy 
an J-A/2 J -a/2 

+ ( X+5X p ■ (0, A t /2) ■ (-AJ2, 1) dx 



+ [ X+6X p ■ (0, -A t /2) ■ (-A x /2, -1) dx 

J X 
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= mA\ x x +Sx + Y(pA t )\y=A/2 + Y(pAt)\v=-A/2 + 0((5x) 2 ) 
« (p-u^ll^ + pAtSx + OdSx) 2 ), (2.8) 

where we have used the smallness of p y to approximate p\ y =±A/2 by p and 
pul by p ■ ul. 

Combining Q-Q, (p]) with: 



jJ n PdV = -(p ■ urA)^ 5 * + 0((5xn (2-9) 
dividing by 5x and sending it to zero, we have: 

(pA) t + (p • u x A) x = 0, 

which is (|1.7| ). 

Next consider i — 1 in the momentum equation, a = x, b = x + 5x. We 
have similarly with ((27 



— / puidV = [ (pu 1 ) t dV+ [ puiJ t dV 
at Jn(t) Jn{t) JQ(t ) 

= — / pu\u-ndS+ \ &i j ■ nj dS + ~pui A t 5x + 0(5x 2 ). (2.10) 

Jdfl(t) JdQ(t) 

We calculate the integrals of ( |2.10| ) below. 
d r 

— / pu x dV = (puiA) t Sx + 0{{5x) 2 ) » (p ■ uA) t ■ 5x + 0((5x) 2 ), (2.11) 
where Mi = on the upper and lower boundaries is used. 

/ pu x u-ndS ={-p-u\A)\ x x +5x + 0(5xp 1 ' 2 ), (2.12) 

Jdtt 

where the smallness of U\ )V in the interior and small width of boundary layer 
0(p}l 2 ) gives the 0(/i 1//2 ) for approximating u\ by ui ■ u\. 



/ -pdijUjdS « -pA\ x x +dx + pA x dx 

Jdtt Jx 

= -pA\ x x +5x +pA x 5x + 0{{5x) 2 ) 
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Noticing that: 

d u = 2p(u 1<x - (u 1>x + u 2 , y )/3), d 12 = 2p(u hy + u 2 , x ). 
It follows that 

-j- 4 _ 2pA t 

Thus the contribution from the left and right boundaries located at x and 
x + Sx is: 

E / dunx = AdT^ = U^ x \t +5x ~ 2j ^\T Sx . (2.13) 

The contribution from the upper and lower boundaries is: 

d u nxdS = -d 11 A x 5x/2\ y=A/2 - d 11 A x 5x/2\ y= _ A/2 

± Jy=±A/2 

= p5xJ20(d y u)\ y=±A/2 . (2.14) 
± 

Similarly, 

J2 di 2 n 2 dS = /j,5xJ20(d y u)\y=±A/2- (2.15) 

± Jy=±A/2 ± 

Since d y u\ y= ±A/2 = 0(/i _1//2 ), the viscous flux from the boundary layers 
are 0(5xp 1 ^ 2 ), much larger than the averaged viscous term 5 x -^(Aui x ) x = 
0{5xp). We notice that the vertically averaged quantities have little de- 
pendence on the viscous boundary layers unless A is on the order 0(p 1 ^ 2 ). 
Hence the quantities from upper and lower edges in ( 2.14 ) and Q2.15|) , and 
that in ( 2.12 ), should balance themselves. Omitting them altogether, and 
combining remaining terms that involve only ul, ~p in the bulk, we end up 
with (after dividing by 5x and sending it to zero): 

4/i 

(p-WA) t +{p-ufA) x = -(p A) x +A x p+pW A t +-^(Au^ x ) x -2fiA tx /3. (2.16) 

Simplifying fl2.16|) with the continuity equation (|1.7f) , we find that: 

_ , ,_ , A t uj 4p x 2jiA- x A tx 

v>u + ui u\ x — —p x / p H — + — A {Au lx ) x — , (2.17) 

A 3p 3p 

which is (II 
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3 Linear Stability Analysis near a Flat Fold 



In this section, we discuss the existence of a neutral oscillation mode of the 
linearized system around a flat fold in the limit p — > 0. Such a mode provides 
an onset condition of oscillation, see similar analysis for the Titze model in 
[^| and fll4l . It turns out that such a mode exists under a theshold condition 
for system (|1.7| j-( T7TT|) because of the A t u/A term. In our case, it is essential 
that the derivation originated with the no slip boundary condition on the 
fold, which allows enough energy of the background flow to transfer to the 
fold to offset the loss there. This energy transfer mechanism is expounded in 
Titze [^] in his body-cover model. We show with a stability analysis that 
our model system supports such a physical picture of flow induced oscillation. 

We assume that the cross section is rectangular. The system (|1.7|)-( |T~TT| ) 
admits constant steady states: (uo,po, A ), p = K,po , satisfying: 

- (3A + 4w Po = f m , (3.1) 

where f m is a constant so that Aq matches the height of the connecting 
vocal tract. We are interested in conditions leading to the small amplitude 
oscillations near the constant steady states. This is similar to Titze [p2]1 , 
where a lumped ODE is proposed and analyzed for the fold center using 
mucosal wave approximation. However here, we calculate directly from ( |1.7| 
)- (|l.ll| ), and do not make any further modeling assumptions. Since we have 
zero background pressure gradient, our results do not compare directly with 
pifl , though qualitative features remain. Another comparison with is 
performed in the next section under quasi-steady approximation where small 
pressure gradient is present. 

Letting u = u + u, p = p + p, p = p + p, A = A + A, and linearizing 
Q-dTTID with p = 0, we get: 

(A p + Ap )t + (p A u + p u A + u A p) x = 0, (3.2) 

u t + u u x H p x = -£-A t , (3.3) 

Po A 
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A a = ok xx - a A t -I3A + Awp. (3.4) 
Equations (|3.2|) - (|3.3|) are written as: 



A pt + PoA t + PqA u x + p u A x + u A p x = 0, (3.5) 

u t + u Q u x + — = -£-A t . (3.6) 
Po 

Applying the operator d t + Uq8 x on ( |3.5| ) and using (|3.6j), we find: 

(d t + u d x )(A p t + p A t + p u A x + u A p x ) - A p xx + u p A t = 0. (3.7) 

Differentiating ( |1 . 1 1| ) gives: pt = K~fp7~ l pt, or p t = kjPq 1 pt upon lineariz- 
ing at p = po- Similarly, p x = kjPq 1 p x . With these relations, equation (|3.7|) 
becomes: 



T(d t + u d x ) 2 p - p xx + ^^-A 



A 

+ A 1 p (d t + u Q d x ) 2 A = 0, (3.8) 



where: 

1 1 



7-1 «2' 
KlPo C 

with c being the speed of sound at air density po- 

Applying the operator A T(d t + u d x ) 2 — A d xx to both sides of (|3.4|) , we 

get: 

A [Td u + 2Fu d xt + (Pu 2 - 1)9^,] ■ (A tt - ok xx + aA t + f3A) 

= -Awpo[{d t + u d x ) 2 A + u A t }. (3.9) 

Substituting the mode A = A m e im ' x+Xt , m! = n ^-, we end up with the 
following algebraic equation of degree four for A: 

[rA 2 + 2Tu m'Xt + (1 - Tu 2 )m' 2 } ■ [A 2 + am 2 + aX + (3] 

= -AwpoA^lX 2 + 2u im'X - u\w! 2 \ - Awu °P° X , ( 3 .i ) 

Aq 
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or: 

TA 4 + (aT + 2Tu m'i)\ 3 + 
(r(/3 + am' 2 ) + 2aTu m'i + AwpoA^ 1 + (1 - Tu 2 )m' 2 )X 2 
+ (2Tu m'i(f3 + am' 2 ) - am' 2 (Tu 2 - 1) 
+Awp Ao \2u tm') + ^E±)\ 
+ [(0 + am' 2 )(l - Yul)m' 2 + Awp A 1 (-u 2 m' 2 )] = 0. (3.11) 

Proposition 3.1 If 

p ul> a(p u + Tp^^), (3.12) 

( \3. 1 1\ ) has a pair of pure imaginary solution A = ±ir], f] ^ being real, which 
implies the existence of a pair of oscillatory modes to the linearized system 
( \3.Q )-(3.4 ) of the form e ±O m7ra /- L +«7*) ; j or rea \ an( i nonzero numbers m and 



rj. The transition into oscillation is a Hopf bifurcation. 



Proof: Let A = irj in (|3.11 ), where r] is real. The real and imaginary parts 
give respectively: 

IV + 2Tu m'r ] 3 ~[T(am' 2 + f3) + m' 2 (l-Tu 2 ) + ^^]r ] 2 

Aq 

r or '(a i a\ 8wp u m' 

+ [-2Tu m (fj + am ) \r] 

Ao 

+ [m /2 (l - Yul){am' 2 + 0) - *&p!£l ] = 0> (3.13) 



and: 



aT V 3 - 2aTu m' V 2 + m' 2 a(l - Tu 2 ) V + AwU f^ = . (3.14) 



For rj 0, a ^ 0, we have from ( 3.14 ): 



IV + 2Tu m'r ] - m' 2 (l - Tu 2 ) - = 0, 

A a 
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so: 



T] = —u m' ± c\Jm' 2 + AwuqPq/ (A a). (3.15) 
Now we regard the left hand side of (|3.13| ) as a continuous function of m', 



call it F(m'). For \m'\ 3> 1, r\ ~ (— uq ± c)m', direct calculation shows: 



F(m') —m' 2 < 0. 



While for Im'l <C 1, 



and: 



provided: 



F{m) ~TA^{-A^- Tp -^) >0 ' 

^p 0>a f*m +r/3 \ (3-16) 



holds, which is just ( |3.12| ). Under ( |3.12| ), F(m) = has a nonzero real solu- 
tion, hence an oscillatory mode solution exists to (|3.2|)-(|3T^). Finally, notic- 



ing that equations ( |3.13|) - (|3.14]) are invariant under the symmetry transform 
(77, m') — > (—77, ~ m ')i we conclude that the oscillatory modes exist as a pair, 
and oscillation appears as a Hopf bifurcation. 

Remark 3.1 Condition says that the fluid energy must be large enough 

to overcome the fold damping due to a. Without the lower order term A t u/ A, 
the same calculation would show that F(m') = —4wpom' 2 /(TAo) < for all 
m' , implying non-existence of oscillatory mode. Condition ( \3. 1 (\ ) is similar to 
the threshold pressure in Titze 's model fffi^ in the sense that minimum energy 
( analogous to minimum lung pressure ) is proportional to the fold damping co- 
efficient and the prephonatory half width A . 
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Remark 3.2 There is another neutral yet nonoscillatory mode from (jS.lS ) 



( \3.14\ ), namely, r] = 0, and: 



*KI = (1 _ m2)Aq - P, (3.17) 

provided the right hand side expression is positive. As we turn on p > but 
small, the oscillation mode in Proposition 2. 1 will generically be slightly per- 
turbed and yet preserves its oscillatory nature. The calculations are tediuous 
and not shown here. 

4 Quasi-steady Approximation and Relations 
with the Titze Model 

The glottal flow is often regarded as nearly quasisteady and inviscid in the 
bulk |jT7|| , and the Bernoulli's law is adopted as an approximation, p2 



among others. In this approximation, the temporal variation of flow variables 
is considered much slower than that of the fold motion. Now let us drop the 
time derivatives, and viscous term in the flow equations to get: 

( P uA) x = 0, (4.1) 
Px . A t u 

uu x = 1 j-, (4.2) 

while keeping the fold dynamic equation and the equation of state the same. 

Our goal is to derive a closed equation for the cross section area A. Inte- 
grating ( |4.1| ) and ( |4.2| ) in x using the equation of state shows: 

puA = Q , (4.3) 
uV2 + ^=r^ dx+ P , (4.4) 



7-1 J-L A 

where Pq and Qq are constants determined by the flow conditions at the inlet 

Atu 
A ' 



x = —L. Note that without the integral term with — *p, (|4.4|) becomes the 



standard Bernoulli's law. 
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Substituting ( |4.3| ) into ( |4.4j ) gives: 



x A t 



L p A 2 



dx + P . 



(4.5) 



We would obtain a closed equation on A if ( |4.5| ) could be solved for p in 
terms of A. 

Suppose that A undergoes small vibration about a constant state A 
which satisfies: 

Ql 



7-1 

I 7Po _ p 

H T — -M)j 



2pgA2 ' 7 -l 

then (|4.5| ) can be solved for p by perturbation. Letting A = A +A, p = po+p, 
we find: 

Ql . , ~ /* A 



Po ^o" Po^o 



i + 7 pr 2 p - Qo 



A 2 p c 



0. 



or: 



P • 7Po 2 



Ql_ 

.3 42 



Qo i Qo 



A + 



If 7Po > or: 



Po At) A 2 oPo A^poJ-L 



7+1 Qo 

Po > 



.4, 



7^' 



which means large pressure for given Qq and Aq, then: 

- Ql i 



p 



Po ^0 



-L p A 



(4.6) 



(4.7) 



(4.8) 



Now substituting ( |4.8| ) into the fold dynamic A equation, we obtain: 



A u - aA xx + a A t + (3A = c aero [ A t + -^i] , 

-l p A 



where: 



?i- 2 - PoAl 



> 0. 



(4.9) 



(4.10) 



7Po 
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Equation ( |4.9| ) can be written as: 

Att - ok xx + (aA - c aero j_ A) t + {(3- c aero Q Po 1 A Q x )i = 0. (4.11) 
Equation ( |4.11| ) is a linear wave equation with damping and pumping. This 



is easy to see from the energy identity. Assuming that A x (±L,t) = 0, we 
have: 

d f L 



dt 



J_ L dx {^Al + l -A] + {(3- c aero Q oP ^A^)A 2 /2 
-a [ L dx A 2 t + c aero ( [ L dxA t ) 2 /2. (4.12) 

J-L J-L 

We shall require that: 

ft - CaeroQoPo^Q 1 > 0, (4.13) 



which is true if po is sufficiently large for fixed ft, A , and Qq. Then ( [4.1 2|) says 



that the rate of change of the total energy in time depends on the balance 
of the natural fold damping (the negative term with prefactor a) and the 
energy input from the aerodynamic flow (the positive term with prefactor 

We see from (|4. 12| ) that spatially sinusoidal modes like e im ' x will render 



J^ L A t dx = 0, hence they do not sustain lossless temporall oscillations. How- 
ever, there are lossless oscillatory solutions with monotone spatial profiles. 
To show a simple solution of this kind, extend the incoming flow uniformly 
to — oo and modify the energy transfer term to /f^ ^ in Q4.11Q . Let us 
keep the exit location still at x = L. Then seek (p = e ux ip (t), v > 0, (p x = A, 
if)(t) satisfies: 

Vlpu + (aV - C aero )lp t + V{ft - CaeroQopQ 1 Aq 1 - aU 2 )lp = 0. 

Since v~ x is the length scale of spatial decay, v~ l is of order O(L). The 
critical condition for temporal oscillation is av = c aero , to remove the damp- 
ing, or 

a 

Caero — ~j~ i 
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or: 



^ 9 aAl 
4wkQ% > -jA 



The temporal oscillation frequency is y /3 — c aer oQoPo ~~ av2 i which 
is positive for large po and small a. Notice that the fold positions at ±L : 
A + h>e~ uL ip(t), A + ve uL ip{t), are not equal, they are above and below A 
simultaneously. When they are above A , the fold is divergent (positive slope 
along the flow direction); and when they are below A Q , the fold is convergent 
(negative slope along the flow direction). However, the two end points are 
never completely out of phase, i.e, when one is above and the other is below 
A . 

The oscillatory solutions are of very different forms with or without mak- 
ing the quasi-steady approximation, which may underestimate the amount 
of the energy transfer from air flow onto the vocal fold. Nevertheless, the 
above spatially monotone solutions under the quasi-steady approximation 
agree qualitatively with those of the Titze model |22[] and shows explicitly 
the effect of the energy transfer term (fZ L jr)- This term is nonlocal and 
different from the local term in Titze model ||22|| , yet is of the same origin 
and reflects the changing directions of the effective air driving force on open 
and closing cycles, a key observation of the Titze model. 

Both our model and the Titze model capture the positive energy feedback 
from air flow into the fold motion, though they are not of the same form. 
Our model has the advantage of a systematic treatment and of characterizing 
the entire fold shape in the continuum. 
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